#Getting the local path to R libraries
if(Sys.info()[["sysname"]]=="Windows"){
	library(Biobase)
}else{
	locFile =file("../Misc/localization_files.txt","r")
	for(line in readLines(locFile)){
		if(substr(line,0,9) == "Rlib_path"){
			Rlib_path = gsub(" ","",sub("=","",sub("Rlib_path","",line)))
		}
	}
	library(Biobase, lib.loc=Rlib_path)
}



expressionXFile="../TXT/exp1_genevsgene.txt"
expressionYFile="../TXT/exp2_genevsgene.txt"
annotationFile="../TXT/annotation_genevsgene.txt"
tissue_nameFile="../TXT/tissuetype_genevsgene.txt"


#################################################
#REFORMING FROM PHP
#################################################
#loading data
require(Biobase) 
expressionX<-read.table(expressionXFile,sep="\t",row.names=1,header=TRUE)
colnames(expressionX)<-sub("^X","",colnames(expressionX)) 

expressionY<-read.table(expressionYFile,sep="\t",row.names=1,header=TRUE)
colnames(expressionY)<-sub("^X","",colnames(expressionY)) 


tissue_name<-read.table(tissue_nameFile,stringsAsFactors=FALSE)[,1]

annotation<-read.table(annotationFile,sep="\t",row.names=1,header=TRUE)

#checking X set and Y set contain data from same samples and that overlapping probesets have same values
if(!all(rownames(expressionX) == rownames(expressionY)))stop("The sampleNames from expressionX and expressionY were not equal") 
for(XandYProbeset in intersect(colnames(expressionX),colnames(expressionY))){
	print(paste("XandYProbeset: ",XandYProbeset))
	if(!all(expressionX[,XandYProbeset] == expressionY[,XandYProbeset]))stop("probesets found in both expressionX and expressionY have differing expression")
	
	probesetsToKeep<-colnames(expressionY)[!colnames(expressionY)%in%XandYProbeset]
	if(length(probesetsToKeep)>1){
		expressionY[,probesetsToKeep]
	}else{
		expressionY<-data.frame(row.names=rownames(expressionY),expressionY[,probesetsToKeep])
		colnames(expressionY)<-probesetsToKeep
	}
}



expression<-cbind(expressionX, expressionY)
y_genes<-colnames(expressionY)


#forming ExpressionSet
expressionset<- new(
		"ExpressionSet", 
		exprs = t(as.matrix(expression)), 
		experimentData = new("MIAME",title=tissue_name))



#Doing statistics
correlations<-data.frame(pvalue=vector(),correlation=vector())
for(x_gene in colnames(expressionX)){
	for(y_gene in colnames(expressionY)){
		cor_result<-cor.test(exprs(expressionset)[x_gene,],exprs(expressionset)[y_gene,])
		correlations[paste(x_gene,y_gene),]<-c(cor_result[[3]],cor_result[[4]])
	}
}
correlations<-correlations[order(abs(correlations[,"correlation"]),decreasing=TRUE),]




#CREATING PNG
if(nrow(correlations)> 15){
	correlationsForPNG<-correlations[1:15,]
	numberOfRows <- 1 + ceiling(nrow(correlationsForPNG) / 3)
}else{
	correlationsForPNG<-correlations
	numberOfRows <- ceiling(nrow(correlationsForPNG) / 3)
}

png(filename = "../Images/genes_vs_genes.png", width = 1440, height = 480*numberOfRows)
par(omi=c(0.2,0.2,0.2,0.2),mai=c(1,1,1,1))
layout(matrix(1:(3*numberOfRows),ncol=3,byrow=TRUE) )
for(i in 1:nrow(correlationsForPNG)){
	rowEntry<-rownames(correlationsForPNG)[i]
	x_gene<-sub(" .+$","",rowEntry)
	y_gene<-sub("^.+ ","",rowEntry)

	plot.default(
			type="p",
			y=exprs(expressionset[y_gene,]),
			x=exprs(expressionset[x_gene,]),
			ylab=list(paste(y_gene,": ",annotation[y_gene,"genesymbol"],sep=""),cex=1.5),
			xlab=list(paste(x_gene,": ",annotation[x_gene,"genesymbol"],sep=""),cex=1.5),
			main=list(paste("Expression plot of",x_gene,"varying with",y_gene),cex=1.5),
			sub=list(paste("correlation: ",signif(correlationsForPNG[rowEntry,"correlation"],3),"pvalue",signif(correlationsForPNG[rowEntry,"pvalue"],3),sep=" "),cex=1.5)
	)
	mtext(paste("       on Y-axis: ",annotation[y_gene,"genesymbol"]," - ",substr(annotation[y_gene,"genename"],1,40),sep=""),adj=0,padj=-1.4,cex=1.2)
	mtext(paste("       on X-axis: ",annotation[x_gene,"genesymbol"]," - ",substr(annotation[x_gene,"genename"],1,40),sep=""),adj=0,padj=-0.2,cex=1.2)
}
if(nrow(correlations)> 15){
	frame()
	frame()
	mtext("more plots available in pdf file",cex=2)
}
garbage <- dev.off()



#CREATING PDF
pdf("../PDF/genes_vs_genes.pdf", width = 3*5, height = 5*5)
layout(matrix(1:15,ncol=3,byrow=TRUE) )
for(i in 1:nrow(correlations)){
	rowEntry<-rownames(correlations)[i]
	x_gene<-sub(" .+$","",rowEntry)
	y_gene<-sub("^.+ ","",rowEntry)
	plot.default(
			type="p",
			y=exprs(expressionset[y_gene,]),
			x=exprs(expressionset[x_gene,]),
			ylab=paste(y_gene,": ",annotation[y_gene,"genesymbol"],sep=""),
			xlab=paste(x_gene,": ",annotation[x_gene,"genesymbol"],sep=""),
			main=list(paste("Expression plot of",x_gene,"varying with",y_gene),cex=0.7),
			sub=list(paste("correlation: ",signif(correlations[rowEntry,"correlation"],3),"pvalue",signif(correlations[rowEntry,"pvalue"],3),sep=" "),cex=0.7)
	)
	mtext(paste("       on Y-axis: ",annotation[y_gene,"genesymbol"]," - ",annotation[y_gene,"genename"],sep=""),adj=0,padj=-1.4,cex=0.7)
	mtext(paste("       on X-axis: ",annotation[x_gene,"genesymbol"]," - ",annotation[x_gene,"genename"],sep=""),adj=0,padj=-0.2,cex=0.7)
}
garbage <- dev.off()




#CREATING TXT
rawData<-cbind(annotation[featureNames(expressionset),],exprs(expressionset))
write.table(rawData, file="../TXT/raw_data.xls", sep="\t",col.names=NA)
